import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from scipy.stats import zscore
import chinese_calendar as calendar
from chinese_calendar import is_workday
#step2 数据分析
# ----------------------------
# 数据加载与列名生成（核心修复）
# ----------------------------

# 加载负荷数据
data = pd.read_csv('data/Area1_Load_clean.csv')
print("原始数据前5行：")
print(data.head())

# 生成正确的负荷列名（T0000~T2345，每15分钟一个，共96列）
def get_load_columns():
    cols = []
    for hour in range(24):  # 小时：00~23
        for minute in [0, 15, 30, 45]:  # 分钟：每15分钟一个间隔
            # 格式化为"小时+分钟"（如00:00→T0000，00:15→T0015，01:00→T0100）
            col_name = f'T{hour:02d}{minute:02d}'
            cols.append(col_name)
    return cols[:96]  # 确保只取96列

# 获取负荷列名
load_cols = get_load_columns()
# 验证列名是否存在于数据中（调试用）
missing_cols = [col for col in load_cols if col not in data.columns]
if missing_cols:
    print(f"警告：数据中缺少以下列名：{missing_cols}")
else:
    print("所有负荷列名均存在于数据中")

cleaned_data= data

# ----------------------------
# 节假日属性处理
# ----------------------------

# 对日期进行处理
cleaned_data['YMD'] = pd.to_datetime(cleaned_data['YMD'])
cleaned_data['weekday'] = cleaned_data['YMD'].dt.weekday
cleaned_data['is_workday'] = cleaned_data['YMD'].apply(is_workday).astype(int)

# 节假日权重设置
Holiday = {
    "New Year's Day": 10,
    "Spring Festival": 30,
    "Tomb-sweeping Day": 10,
    "Labour Day": 10,
    "Dragon Boat Festival": 10,
    "National Day": 20,
    "Mid-autumn Festival": 10
}

Holiday_Info = []
Holiday_Name = []
time = pd.to_datetime(cleaned_data['YMD'])

for t_date in time.dt.date:
    on_holiday, holiday_name = calendar.get_holiday_detail(t_date)
    if on_holiday:
        # 优先使用自定义权重，周末默认5
        t_holiday = Holiday.get(holiday_name, 5 if t_date.weekday() >= 5 else 0)
        Holiday_Name.append(holiday_name)
    else:
        t_holiday = 0
        Holiday_Name.append('None')
    Holiday_Info.append(t_holiday)

# 添加节假日属性列
if 'Holiday_Info' in cleaned_data.columns:
    cleaned_data.drop('Holiday_Info', axis=1, inplace=True)
if 'Holiday_Name' in cleaned_data.columns:
    cleaned_data.drop('Holiday_Name', axis=1, inplace=True)

cleaned_data.insert(102, 'Holiday_Info', Holiday_Info)
cleaned_data.insert(103, 'Holiday_Name', Holiday_Name)

# 保存含节假日属性的数据
cleaned_data.to_csv('data/Area1_Load_clean.csv', index=False)
print("含节假日属性的负荷数据已更新保存")

# ----------------------------
# 季节性负荷可视化
# ----------------------------

# 选取2011年春夏秋冬各一天的数据
data = pd.read_csv('data/Area1_Load_clean.csv')
data['YMD'] = pd.to_datetime(data['YMD'])

spring = data[data['YMD'] == '2011-04-15']
summer = data[data['YMD'] == '2011-07-15']
autumn = data[data['YMD'] == '2011-10-15']
winter = data[data['YMD'] == '2012-01-16']

# 确保数据存在
if len(spring) == 0:
    print("警告：未找到2011-04-15的数据，使用首条春季数据替代")
    spring = data[(data['YMD'] >= '2011-03-01') & (data['YMD'] <= '2011-05-31')].head(1)

plt.figure(figsize=(12, 6))
plt.plot(range(len(spring.iloc[0, 1:97])), spring.iloc[0, 1:97], linestyle='-', label='Spring')
plt.plot(range(len(summer.iloc[0, 1:97])), summer.iloc[0, 1:97], linestyle='-', label='Summer')
plt.plot(range(len(autumn.iloc[0, 1:97])), autumn.iloc[0, 1:97], linestyle='-', label='Autumn')
plt.plot(range(len(winter.iloc[0, 1:97])), winter.iloc[0, 1:97], linestyle='-', label='Winter')
plt.xlabel('Time/15min')
plt.ylabel('Load/MW')
plt.legend()
plt.title('Daily Load in Different Seasons (2011)')
plt.show()

# ----------------------------
# 每周负荷趋势可视化
# ----------------------------

# 选取2011年各季节的完整一周（周一至周日）
spring_start = pd.to_datetime('2011-04-11')  # 周一
summer_start = pd.to_datetime('2011-07-11')
autumn_start = pd.to_datetime('2011-10-10')
winter_start = pd.to_datetime('2012-01-09')

# 验证起始日期为周一
assert spring_start.weekday() == 0, "春季起始日期不是周一"
assert summer_start.weekday() == 0, "夏季起始日期不是周一"
assert autumn_start.weekday() == 0, "秋季起始日期不是周一"
assert winter_start.weekday() == 0, "冬季起始日期不是周一"

# 提取每周数据
spring_week = data[(data['YMD'] >= spring_start) & (data['YMD'] < spring_start + pd.Timedelta(days=7))]
summer_week = data[(data['YMD'] >= summer_start) & (data['YMD'] < summer_start + pd.Timedelta(days=7))]
autumn_week = data[(data['YMD'] >= autumn_start) & (data['YMD'] < autumn_start + pd.Timedelta(days=7))]
winter_week = data[(data['YMD'] >= winter_start) & (data['YMD'] < winter_start + pd.Timedelta(days=7))]

# 展平为连续时间序列
spring_week = spring_week.iloc[:, 1:97].values.flatten()
summer_week = summer_week.iloc[:, 1:97].values.flatten()
autumn_week = autumn_week.iloc[:, 1:97].values.flatten()
winter_week = winter_week.iloc[:, 1:97].values.flatten()

# 可视化每周负荷趋势
plt.figure(figsize=(12, 6))
plt.plot(range(len(spring_week)), spring_week, linestyle='-', label='Spring Week')
plt.plot(range(len(summer_week)), summer_week, linestyle='-', label='Summer Week')
plt.plot(range(len(autumn_week)), autumn_week, linestyle='-', label='Autumn Week')
plt.plot(range(len(winter_week)), winter_week, linestyle='-', label='Winter Week')
plt.xlabel('Time/15min')
plt.ylabel('Load/MW')
plt.legend()
plt.title('Weekly Load Trends in Different Seasons')
plt.show()

# ----------------------------
# 工作日/周末/节假日负荷对比
# ----------------------------

# 选取样本日期
workday = data[data['YMD'] == '2011-04-11']  # 周一（工作日）
weekend = data[data['YMD'] == '2011-04-16']   # 周六（周末）
holiday = data[(data['Holiday_Name'] == 'Spring Festival') & (data['YMD'] == '2011-02-02')]

# 可视化对比
plt.figure(figsize=(12, 6))
plt.plot(range(len(workday.iloc[0, 1:97])), workday.iloc[0, 1:97], linestyle='-', label='Workday')
plt.plot(range(len(weekend.iloc[0, 1:97])), weekend.iloc[0, 1:97], linestyle='-', label='Weekend')
# plt.plot(range(len(holiday.iloc[0, 1:97])), holiday.iloc[0, 1:97], linestyle='-', label='Spring Festival')
plt.xlabel('Time/15min')
plt.ylabel('Load/MW')
plt.legend()
plt.title('Load Comparison: Workday vs Weekend')
plt.show()

# ----------------------------
# 温度与负荷相关性分析（修改后）
# ----------------------------

# 加载数据并筛选时间范围
data_temp = pd.read_csv('data/Area1_Load_clean.csv')
weather_temp = pd.read_csv('data/Area1_Weather_clean.csv')

# 转换日期格式
data_temp['YMD'] = pd.to_datetime(data_temp['YMD'])
weather_temp['YMD'] = pd.to_datetime(weather_temp['YMD'])

# 关键修复：通过日期合并数据，确保两者日期完全对齐
# 只保留双方都存在的日期
merged_data = pd.merge(
    data_temp[['YMD', 'MeanLoad']],  # 负荷数据只保留日期和平均负荷
    weather_temp[['YMD', 'Max_Temperature', 'Min_Temperature', 'Avg_Temperature']],  # 气象数据保留日期和温度
    on='YMD',  # 按日期合并
    how='inner'  # 内连接：只保留双方都有的日期
)

# 筛选2012-2014年数据（在合并后筛选，确保日期一致）
time_mask = (merged_data['YMD'] >= '2012-01-01') & (merged_data['YMD'] < '2015-01-01')
merged_data = merged_data[time_mask]

# 处理温度异常值（>50°C用前一天替换）
for col in ['Max_Temperature', 'Avg_Temperature', 'Min_Temperature']:
    merged_data[col] = np.where(
        merged_data[col] > 50,
        merged_data[col].shift(1),  # 用前一天数据替换
        merged_data[col]
    )

# 可视化温度与负荷关系（此时x和y长度一定相同）
plt.figure(figsize=(12, 6))
plt.scatter(merged_data['Avg_Temperature'], merged_data['MeanLoad'], marker='*', label='Avg Temperature')
plt.scatter(merged_data['Max_Temperature'], merged_data['MeanLoad'], marker='o', label='Max Temperature', alpha=0.5)
plt.scatter(merged_data['Min_Temperature'], merged_data['MeanLoad'], marker='s', label='Min Temperature', alpha=0.5)
plt.xlabel('Temperature/°C')
plt.ylabel('Daily Mean Load/MW')
plt.legend()
plt.title('Relationship Between Temperature and Load')
plt.show()

# 保存处理后的气象数据（可选）
weather_temp.to_csv('data/Area1_Weather_clean.csv', index=False)

# ----------------------------
# 降雨量与负荷相关性分析
# ----------------------------

# 筛选2013年4月数据
data_rain = pd.read_csv('data/Area1_Load_clean.csv')
weather_rain = pd.read_csv('data/Area1_Weather_clean.csv')

data_rain['YMD'] = pd.to_datetime(data_rain['YMD'])
weather_rain['YMD'] = pd.to_datetime(weather_rain['YMD'])

april_mask = (weather_rain['YMD'] >= '2013-04-01') & (weather_rain['YMD'] < '2013-05-01')
weather_rain = weather_rain[april_mask]

# 选取不同降雨强度的样本
sunny = weather_rain[weather_rain['Rainfall'] == 0].sample(1) if not weather_rain[weather_rain['Rainfall'] == 0].empty else None
light_rain = weather_rain[(weather_rain['Rainfall'] > 0) & (weather_rain['Rainfall'] <= 10)].sample(1) if not weather_rain[(weather_rain['Rainfall'] > 0) & (weather_rain['Rainfall'] <= 10)].empty else None
moderate_rain = weather_rain[(weather_rain['Rainfall'] > 10) & (weather_rain['Rainfall'] <= 25)].sample(1) if not weather_rain[(weather_rain['Rainfall'] > 10) & (weather_rain['Rainfall'] <= 25)].empty else None
heavy_rain = weather_rain[weather_rain['Rainfall'] > 25].sample(1) if not weather_rain[weather_rain['Rainfall'] > 25].empty else None

# 可视化降雨与负荷关系
plt.figure(figsize=(12, 6))
if sunny is not None:
    sunny_load = data_rain[data_rain['YMD'] == sunny['YMD'].values[0]].iloc[0, 1:97].values
    plt.plot(range(len(sunny_load)), sunny_load, label='Sunny (0mm)')
if light_rain is not None:
    light_load = data_rain[data_rain['YMD'] == light_rain['YMD'].values[0]].iloc[0, 1:97].values
    plt.plot(range(len(light_load)), light_load, label='Light Rain (<=10mm)')
if moderate_rain is not None:
    moderate_load = data_rain[data_rain['YMD'] == moderate_rain['YMD'].values[0]].iloc[0, 1:97].values
    plt.plot(range(len(moderate_load)), moderate_load, label='Moderate Rain (10-25mm)')
if heavy_rain is not None:
    heavy_load = data_rain[data_rain['YMD'] == heavy_rain['YMD'].values[0]].iloc[0, 1:97].values
    plt.plot(range(len(heavy_load)), heavy_load, label='Heavy Rain (>25mm)')

plt.xlabel('Time/15min')
plt.ylabel('Load/MW')
plt.legend()
plt.title('Load Under Different Rainfall Conditions (April 2013)')
plt.show()

# ----------------------------
# 湿度与负荷相关性分析（修改后）
# ----------------------------

data_hum = pd.read_csv('data/Area1_Load_clean.csv')
weather_hum = pd.read_csv('data/Area1_Weather_clean.csv')

data_hum['YMD'] = pd.to_datetime(data_hum['YMD'])
weather_hum['YMD'] = pd.to_datetime(weather_hum['YMD'])

# 按日期合并，确保数据对齐
merged_hum = pd.merge(
    data_hum[['YMD', 'MeanLoad']],
    weather_hum[['YMD', 'Avg_Humidity']],
    on='YMD',
    how='inner'
)

# 筛选2012-2014年数据
hum_mask = (merged_hum['YMD'] >= '2012-01-01') & (merged_hum['YMD'] < '2015-01-01')
merged_hum = merged_hum[hum_mask]

# 可视化湿度与负荷关系
plt.figure(figsize=(12, 6))
plt.scatter(merged_hum['Avg_Humidity'], merged_hum['MeanLoad'], marker='*', label='Relative Humidity')
plt.xlabel('Relative Humidity/%')
plt.ylabel('Daily Mean Load/MW')
plt.legend()
plt.title('Relationship Between Humidity and Load')
plt.show()

# ----------------------------
# ----------------------------
# 气象因素与负荷相关性热力图（修改后）
# ----------------------------

# 加载并合并数据（确保日期对齐）
data_corr = pd.read_csv('data/Area1_Load_clean.csv')
weather_corr = pd.read_csv('data/Area1_Weather_clean.csv')

# 转换日期格式
data_corr['YMD'] = pd.to_datetime(data_corr['YMD'])
weather_corr['YMD'] = pd.to_datetime(weather_corr['YMD'])

# 按日期合并数据（内连接，只保留双方都有的日期）
merged_corr = pd.merge(
    data_corr[['YMD', 'MeanLoad']],  # 负荷数据：日期+平均负荷
    weather_corr[['YMD', 'Max_Temperature', 'Min_Temperature', 'Avg_Temperature', 'Avg_Humidity', 'Rainfall']],  # 气象数据
    on='YMD',
    how='inner'
)

# 筛选2012-2014年数据
time_mask = (merged_corr['YMD'] >= '2012-01-01') & (merged_corr['YMD'] < '2015-01-01')
merged_corr = merged_corr[time_mask]

# 计算相关系数
correlation = merged_corr[['Max_Temperature', 'Min_Temperature', 'Avg_Temperature', 'Avg_Humidity', 'Rainfall', 'MeanLoad']].corr()

# 绘制热力图
plt.figure(figsize=(10, 6))
sns.heatmap(correlation, annot=True, cmap='coolwarm', fmt='.3f', square=True)
plt.title('Correlation Between Meteorological Factors and Load')
plt.show()

# ----------------------------
# 负荷数据小时级聚合（修改后）
# ----------------------------

data_hour = pd.DataFrame()
data = pd.read_csv('data/Area1_Load_clean.csv')
data['YMD'] = pd.to_datetime(data['YMD'])

for i in range(len(data)):
    t = data['YMD'][i]
    # 每4列取1列（15min→1h）
    load = data.iloc[i, 1:97:4].values
    temp = pd.DataFrame({
        'Time': pd.date_range(t, periods=24, freq='h'),
        'Load': load
    })
    data_hour = pd.concat([data_hour, temp], ignore_index=True)

data_hour.to_csv('data/Area1_Load_hour.csv', index=False)
print("小时级负荷数据已保存至 data/Area1_Load_hour.csv")

# ----------------------------
# 负荷自相关分析
# ----------------------------

data_acf = pd.read_csv('data/Area1_Load_hour.csv')
data_acf['Time'] = pd.to_datetime(data_acf['Time'])

# 计算自相关系数（滞后1-240小时）
autocorrelation = [data_acf['Load'].autocorr(lag=i) for i in range(1, 241)]

# 可视化自相关
plt.figure(figsize=(12, 6))
plt.plot(range(1, 241), autocorrelation, marker='*', linestyle='-')
plt.xlabel('Lag/hour')
plt.ylabel('Autocorrelation Coefficient')
plt.title('Autocorrelation of Hourly Load')
plt.grid(alpha=0.3)
plt.show()

# ----------------------------
# 负荷数据标准化处理
# ----------------------------

# Min-Max标准化
data_minmax = pd.read_csv('data/Area1_Load_clean.csv')
load_cols = get_load_columns()  # 重新获取列名

# 对负荷列标准化
for col in load_cols + ['MaxLoad', 'MinLoad', 'MeanLoad']:
    data_minmax[col] = (data_minmax[col] - data_minmax[col].min()) / (data_minmax[col].max() - data_minmax[col].min() + 1e-6)

data_minmax.to_csv('data/Area1_Load_minmax.csv', index=False)
print("Min-Max标准化数据已保存至 data/Area1_Load_minmax.csv")

# Z-Score标准化
data_zscore = pd.read_csv('data/Area1_Load_clean.csv')
for col in load_cols + ['MaxLoad', 'MinLoad', 'MeanLoad']:
    mean = data_zscore[col].mean()
    std = data_zscore[col].std()
    data_zscore[col] = (data_zscore[col] - mean) / (std + 1e-6)  # 避免除以零

data_zscore.to_csv('data/Area1_Load_zscore.csv', index=False)
print("Z-Score标准化数据已保存至 data/Area1_Load_zscore.csv")